home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / DBio / DSeqList.cpp < prev    next >
Encoding:
Text File  |  1995-12-17  |  13.3 KB  |  596 lines  |  [TEXT/R*ch]
  1. // DSeqList.cp
  2. // d.g.gilbert
  3.  
  4.  
  5. #include "DSeqList.h"
  6. #include "DSequence.h"
  7. #include "DSeqFile.h"
  8. #include <ncbi.h>
  9. #include <dgg.h>
  10. #include <DList.h>
  11. #include <DFile.h>
  12. #include <ureadseq.h>
  13.  
  14.  
  15.  
  16. Global short gLinesPerSeqWritten = 0;
  17. short DSeqList::gMinCommonPercent = 70;
  18.  
  19. // DSeqList  -------------------------------------
  20.  
  21. DSeqList::DSeqList() :
  22.     DList(),
  23.     fSortOrder(kSortByItem)
  24. {
  25. }
  26.  
  27. DSeqList::~DSeqList()
  28. {
  29.         // PLUG MEMORY LEAK !! ~DList doesn't call object destructors !!!!
  30.     if (fDeleteObjects) {
  31.         long i, n= GetSize();
  32.         for (i=0; i<n; i++) {
  33.             DSequence* aseq= this->SeqAt(i);
  34.             delete aseq;
  35.             }
  36.         fDeleteObjects= false;
  37.         }
  38. }
  39.  
  40. Boolean DSeqList::IsEmpty()
  41. {
  42.     if (GetSize() > 1) 
  43.         return FALSE;
  44.     else if (GetSize() == 1) {
  45.         DSequence* aSeq= (DSequence*) this->First();
  46.         return (aSeq && aSeq->LengthF() > 0);
  47.         }    
  48.     else
  49.       return TRUE;
  50. }
  51.  
  52.  
  53.  
  54.  
  55.  
  56. class CFileIndex 
  57. {
  58. public:
  59.     long*    fIndices;
  60.     long    fMax, fNum;
  61.  
  62.     CFileIndex();
  63.     ~CFileIndex();
  64.      long* Indices() { return fIndices; }
  65.      long    IndexCount() { return fNum; }
  66.     void Indexit( long index);
  67.     void Indexit( DFile* aFile);
  68.     void IndexEOF( DFile* aFile);
  69. };
  70.  
  71. CFileIndex::CFileIndex() 
  72. {
  73.     fNum= 0;
  74.     fMax= 0;
  75.     fIndices= (long*) MemGet(sizeof(long),true);
  76. }
  77.  
  78. CFileIndex::~CFileIndex()
  79. {
  80.     if (fIndices) MemFree( fIndices);
  81. }
  82.      
  83. void CFileIndex::Indexit( long index)
  84. {
  85.     if (fNum >= fMax){
  86.         fMax += 20; 
  87.         fIndices= (long*) MemMore( fIndices, sizeof(long)*fMax); 
  88.         }
  89.     if (fIndices) fIndices[fNum++]= index;
  90. }
  91.  
  92. void CFileIndex::Indexit( DFile* aFile)
  93. {
  94.     Indexit(aFile->Tell());
  95. }
  96.  
  97. void CFileIndex::IndexEOF( DFile* aFile)
  98. {
  99.     Indexit(aFile->Tell());
  100.     fNum--;
  101. }
  102.     
  103.  
  104.  
  105.     
  106.  
  107. void DSeqList::DoWrite( DFile* aFile, short format)  // change to iostreams !!
  108. {
  109.     Boolean        needSameSize = false, sizesDiffer = false, isInterleaved = false;
  110.     short         iseq, nseqs;
  111.     short            seqtype = DSequence::kOtherSeq;
  112.     long          minbases = -1;
  113.     long             start, nbases;
  114.     DSequence * aSeq;
  115.     DFile         * outfile;
  116.     DFile            * tmpfile = NULL;
  117.     CFileIndex aFileIndex;
  118.  
  119.     
  120.     nseqs= this->GetSize();
  121.     outfile= aFile;
  122.     
  123.     if (nseqs>1) // deal w/ one-per-file formats 
  124.         switch (format) {
  125.             case DSeqFile::kGCG            : format= DSeqFile::kMSF; break;
  126.             case DSeqFile::kStrider    : format= DSeqFile::kIG; break;
  127.             case DSeqFile::kNoformat:
  128.             case DSeqFile::kPlain:
  129.             case DSeqFile::kUnknown : format= DSeqFile::kGenBank; break;
  130.             default                :    break;
  131.             }
  132.         
  133.     //this->Each( findMinSize);    
  134.     for (iseq= 0; iseq<nseqs; iseq++) {
  135.         aSeq= this->SeqAt(iseq);
  136.         aSeq->GetSelection( start, nbases);
  137.         seqtype= aSeq->Kind();
  138.         //if (minbases<0) minbases= nbases;
  139.         if (minbases<=0) minbases= nbases; // ?? temp fix problem w/ empty seq
  140.         else {
  141.             if (nbases!=minbases) sizesDiffer= true;
  142.             minbases= Min( nbases, minbases);
  143.             }
  144.         }
  145.         
  146.     if (!gOutputName) gOutputName= (char*)aFile->GetShortname();
  147.     DSeqFile::WriteSeqHeader( aFile, format, nseqs, minbases, gOutputName, needSameSize, isInterleaved);
  148.     gOutputName= NULL;
  149.     
  150.     if (isInterleaved) {
  151.         tmpfile= new DTempFile();
  152.         outfile= tmpfile;
  153.         }
  154.     
  155.     if ((needSameSize || isInterleaved) && sizesDiffer) {
  156.         //this->Each( SetSameSize);
  157.         for (iseq= 0; iseq<nseqs; iseq++) {
  158.             aSeq= this->SeqAt(iseq);
  159.             aSeq->GetSelection( start, nbases);
  160.             aSeq->SetSelection( start, minbases); 
  161.             }
  162.         }
  163.  
  164.     if (gPretty.isactive && gPretty.numtop) {
  165.         aSeq= this->SeqAt(nseqs-1);
  166.         gPretty.numline = 1;
  167.         if (isInterleaved) aFileIndex.Indexit(outfile);
  168.         aSeq->DoWriteSelection( outfile, format);
  169.         gPretty.numline = 2;
  170.         if (isInterleaved) aFileIndex.Indexit(outfile);
  171.         aSeq->DoWriteSelection( outfile, format);
  172.         gPretty.numline = 0;
  173.         }
  174.  
  175.     for (iseq= 0; iseq<nseqs; iseq++) {
  176.         aSeq= this->SeqAt(iseq);
  177.         if (isInterleaved) aFileIndex.Indexit(outfile);
  178.         aSeq->DoWriteSelection( outfile, format);
  179.         }
  180.                 
  181.     if (gPretty.isactive && gPretty.numbot) {
  182.         aSeq= this->SeqAt(nseqs-1);
  183.         gPretty.numline = 2;
  184.         if (isInterleaved) aFileIndex.Indexit(outfile);
  185.         aSeq->DoWriteSelection( outfile, format);
  186.         gPretty.numline = 1;
  187.         if (isInterleaved) aFileIndex.Indexit(outfile);
  188.         aSeq->DoWriteSelection( outfile, format);
  189.         gPretty.numline = 0;
  190.         }
  191.  
  192.     if (isInterleaved) {
  193.         aFileIndex.IndexEOF( outfile);
  194.         //tmpfile->Close();   // don't need
  195.         //tmpfile->Open("r"); // don't need
  196.         DSeqFile::DoInterleave( aFile, tmpfile, format, seqtype, nseqs,
  197.                                     minbases, gLinesPerSeqWritten, 
  198.                                     aFileIndex.IndexCount(), aFileIndex.Indices());
  199.         outfile= aFile;
  200.         //tmpfile->Close(); // don't need
  201.         tmpfile->Delete();  // debug tmp
  202.         }
  203.         
  204.     DSeqFile::WriteSeqTrailer( aFile, format, this->GetSize(), gLinesPerSeqWritten, minbases);
  205. }
  206.  
  207.  
  208.  
  209. void DSeqList::DoWrite( char* aFileName, short format)  
  210. {
  211.     gOutputName= aFileName;
  212.     DFile* aFile= new DFile( aFileName, "w");
  213.     aFile->Open();
  214.     DoWrite( aFile, format);
  215.     aFile->Close();
  216. }
  217.  
  218.  
  219. #if 0
  220. Handle DSeqList::doWriteHandle( integer format)  
  221. VAR
  222.         Handle        aHand, handSave;
  223.         Boolean        needSameSize, sizesDiffer, isInterleaved;
  224.         integer        aRow, nseqs, seqtype, err, vref, tempFileRef;
  225.         longint        minbases, dirID;
  226.         
  227.     void findMinSize( DSequence* aSeq)
  228.     longint        var  start, nbases;
  229.     {
  230.         aSeq->GetSelection( start, nbases);
  231.         seqtype= aSeq->fKind;
  232.         if (minbases<0) then minbases= nbases
  233.         else {
  234.             if (nbases!=minbases) then sizesDiffer= true;
  235.             minbases= min( nbases, minbases);
  236.             }
  237.     }
  238.     
  239.     void SetSameSize( DSequence* aSeq)
  240.     longint        var  start, nbases;
  241.     {
  242.         aSeq->GetSelection( start, nbases);
  243.         aSeq->SetSelection( start, minbases);
  244.     }
  245.     
  246.     void writeRangeToHand( DSequence* aSeq)
  247.     longint        VAR index;
  248.     {
  249.         if (isInterleaved) then begin
  250.             if (0!=GetFPos( tempFileRef, index)) then ;
  251.             Indexit(index);
  252.             aSeq->doWriteSelection( tempFileRef, format);
  253.             end
  254.         else
  255.             aSeq->doWriteSelectionHandle( aHand, format);
  256.     }
  257.         
  258.     FailInfo        VAR  fi;
  259.     void HndlFailure(OSErr error, long message)        
  260.     {
  261.         IndexCleanup(); 
  262.         aHand= DisposeIfHandle( aHand);
  263.     }
  264.     
  265.     
  266. {
  267.     gOutindex= NULL; //make sure for failure
  268.     aHand= NULL;
  269.     minbases= -1; 
  270.     needSameSize= false;
  271.     isInterleaved= false;
  272.     sizesDiffer= false;
  273.     seqtype= kOtherSeq;
  274.     nseqs= this->GetSize();
  275.     this->Each( findMinSize);        
  276.  
  277.     if (nseqs>1)) // deal w/ one-per-file formats 
  278.         switch (format) {
  279.             kGCG        : format= kMSF;
  280.             kStrider: format= kIG;
  281.             kNoformat, kPlain, kUnknown: format= kGenbank;
  282.             }
  283.  
  284.     CatchFailures(fi, HndlFailure);
  285.     aHand= NewPermHandle(0);
  286.     FailMemError();
  287.     
  288.     DSeqFile::WriteSeqHeader(longint(aHand), TRUE, format, nseqs, minbases, gOutputName,
  289.             needSameSize, isInterleaved);
  290.             
  291.     if (isInterleaved) {
  292.         IndexSetup();
  293.         err= FindFolder(kOnSystemDisk,kTemporaryFolderType,kCreateFolder,vRef,dirID);
  294.         if (err!=0)) { vref= gAppVolRef; dirID= gAppDirID; }
  295.         if 0!=HDelete(vRef, dirID, kInterleaveTemp) then ;
  296.         FailOSErr( HCreate( vRef, dirID, kInterleaveTemp, kSappSig, 'TEXT'));
  297.         FailOSErr( HOpen( vRef, dirID, kInterleaveTemp, fsRdWrPerm, tempFileRef));
  298.         }
  299.     
  300.     if (((needSameSize || isInterleaved/*?*/) && sizesDiffer)) {
  301.         /*---
  302.         ParamText('all sequences must have same # bases for this format', '', 'WriteSeq','');    
  303.         Failure( -1, msgMyError); 
  304.         ---*/
  305.         this->Each( SetSameSize);
  306.         }
  307.  
  308.     this->Each( writeRangeToHand);        
  309.  
  310.     if (isInterleaved) {
  311.         indexEOF( tempFileRef);
  312.         DoInterleave( longint(aHand), TRUE/*isHandle*/, format, seqtype, nseqs,
  313.                                     minbases, gLinesPerSeqWritten, 
  314.                                     gNoutindex, Handle(gOutindex), tempFileRef);
  315.  
  316.         FailOSErr( FSClose(tempFileRef));
  317.         if 0!=HDelete(vRef, dirID, kInterleaveTemp) then ;
  318.         IndexCleanup();
  319.         }
  320.  
  321.     WriteSeqTrailer(longint(aHand), TRUE, format, this->GetSize(), gLinesPerSeqWritten, minbases);
  322.  
  323.     Success(fi);    
  324.     doWriteHandle= aHand;
  325. }
  326.  
  327. #endif
  328.  
  329.  
  330. void DSeqList::AddNewSeq()
  331. {
  332.     DSequence* aSeq = new DSequence();
  333.     InsertLast( aSeq);
  334.     //aSeq->fIndex= this->GetSize();
  335. }
  336.  
  337.  
  338. void DSeqList::ClearSelections()
  339. {
  340.     long i, nseq= GetSize();
  341.     for (i= 0; i<nseq; i++)  
  342.         SeqAt(i)->ClearSelection();
  343. }
  344.  
  345. short DSeqList::ZeroOrigin()
  346. {
  347.     DSequence* aSeq;
  348.     short shiftall= 0;
  349.     long i, nseq= GetSize();
  350.     
  351.     for (i= 0; i<nseq; i++) {
  352.         aSeq= SeqAt(i);
  353.         if (aSeq->Origin()<0) shiftall= Min(shiftall, aSeq->Origin());
  354.         }
  355.         
  356.     for (i= 0; i<nseq; i++) {
  357.         aSeq= SeqAt(i);
  358.         aSeq->SetOrigin(aSeq->Origin() - shiftall); // -= shiftall;
  359.         if (aSeq->Origin() > 0) {
  360.             aSeq->InsertSpacers( 0, aSeq->Origin(), DSequence::indelEdge);
  361.             aSeq->SetOrigin(0);
  362.             }
  363.         }
  364.     return shiftall;
  365. }
  366.  
  367.  
  368.  
  369. short dslCompareByName( void* item1, void* item2)
  370. {
  371.     char*    name1= (char*) ((DSequence*)item1)->Name();
  372.     char*    name2= (char*) ((DSequence*)item2)->Name();
  373.     return StrICmp( name1, name2);
  374. }
  375.  
  376. short dslCompareBySize( void* item1, void* item2)
  377. {
  378.         // inverted sort, largest 1st
  379.     long name1= ((DSequence*)item1)->LengthF();
  380.     long name2= ((DSequence*)item2)->LengthF();
  381.     if (name1 > name2) return -1;
  382.     else if (name1 < name2) return 1;
  383.     else  return 0;
  384. }
  385.  
  386. short dslCompareByKind( void* item1, void* item2)
  387. {
  388.     short name1= ((DSequence*)item1)->Kind();
  389.     short name2= ((DSequence*)item2)->Kind();
  390.     if (name1 > name2) return 1;
  391.     else if (name1 < name2) return -1;
  392.     else  return 0;
  393. }
  394.  
  395. short dslCompareByDate( void* item1, void* item2)
  396. {
  397.         // inverted sort, latest 1st
  398.     unsigned long name1= ((DSequence*)item1)->ModTime();
  399.     unsigned long name2= ((DSequence*)item2)->ModTime();
  400.     if (name1 > name2) return -1;
  401.     else if (name1 < name2) return 1;
  402.     else  return 0;
  403. }
  404.  
  405.  
  406. Boolean DSeqList::SortList(Sorts sortorder)
  407. {
  408.     if (sortorder == fSortOrder)
  409.         return false;
  410.     else {
  411.         fSortOrder= sortorder;
  412.         switch( sortorder) {
  413.             case kSortByKind :
  414.                 this->SortBy( dslCompareByKind);
  415.                 break;
  416.             case kSortBySize : 
  417.                 this->SortBy( dslCompareBySize);
  418.                 break;
  419.             case kSortByDate :
  420.                 this->SortBy( dslCompareByDate);
  421.                 break;
  422.             case kSortByName : 
  423.                 this->SortBy( dslCompareByName);
  424.                 break;
  425.             case kSortByItem: 
  426.             default:
  427.                 this->Sort();
  428.                 break;
  429.             }
  430.         return true;
  431.         }
  432. }
  433.  
  434. DSequence* DSeqList::FindNamedSeq(char* name, Boolean respectCase)
  435. {    
  436.     if (name && *name) {
  437.         DSequence* aSeq;
  438.         long i, nseq= GetSize();
  439.         if (respectCase) for (i= 0; i<nseq; i++) {
  440.             aSeq= SeqAt(i);
  441.           if (Nlm_StringCmp(aSeq->Name(), name)==0) return aSeq;
  442.             }
  443.         else for (i= 0; i<nseq; i++) {
  444.             aSeq= SeqAt(i);
  445.             if (Nlm_StringICmp(aSeq->Name(), name)==0) return aSeq;
  446.             }
  447.         }
  448.     return NULL;    
  449. }
  450.  
  451. DSequence* DSeqList::Consensus(void)
  452. {    
  453.     long i, nseq= GetSize();
  454.     for (i= 0; i<nseq; i++) {
  455.         DSequence* aSeq= SeqAt(i);
  456.         if (aSeq->IsConsensus()) return aSeq;
  457.         }
  458.     return NULL;    
  459. }
  460.  
  461.  
  462. short DSeqList::ConsensusRow()
  463. {
  464.     DSequence* cons= this->Consensus();        
  465.     if (cons) 
  466.         return this->GetIdentityItemNo( cons);
  467.     else
  468.         return -1;
  469. }
  470.  
  471.  
  472.  
  473.     
  474. void DSeqList::MakeConsensus()
  475. {        
  476.     DSequence* cons= this->Consensus();        
  477.     if (!cons) {
  478.         cons = new DSequence(); 
  479.         cons->SetName(DSequence::kConsensus);
  480.         InsertLast( cons); 
  481.         }
  482.         
  483.     if (cons) {
  484.         //short arow= this->GetIdentityItemNo( cons);
  485.         char* hCon= cons->Bases();   
  486.         long conlen= StrLen( hCon); //?? use aSeq->fBaseLength;
  487.         long count= 0;
  488.         long i, iseq, nseq= GetSize();
  489.         for (iseq= 0; iseq<nseq; iseq++) {
  490.             DSequence* aSeq= this->SeqAt(iseq);
  491.             if (!aSeq->IsConsensus() && aSeq->Kind() != DSequence::kOtherSeq) {
  492.                 count++;
  493.                 char* hSeq= aSeq->Bases();
  494.                 long  len = aSeq->LengthF();
  495.                 if (len > conlen) {
  496.                     conlen= len;
  497.                     hCon= (char*) MemMore( hCon, conlen);
  498.                     cons->UpdateBases( hCon, conlen);
  499.                     }
  500.                 if (count==1) {
  501.                     for (i= 0; i<len; i++) hCon[i]= hSeq[i];
  502.                     }            
  503.                 else {
  504.                     for (i= 0; i<len; i++) 
  505. // !! THis is only good for Nucleic bases, not for Aminos !!
  506. // change to table of matches !!
  507.                         hCon[i]= DSequence::NucleicConsensus(
  508.                             DSequence::NucleicBits(hCon[i]), 
  509.                             DSequence::NucleicBits(hSeq[i]));  
  510.                     }
  511.                 }
  512.             }
  513.         }
  514. }
  515.  
  516.  
  517. char* DSeqList::FindCommonBases(short minCommonPerCent, char*& firstCommon)
  518. {  
  519.     enum { kLetrange = 26, kLetlast = 255 };
  520.     char    * hSeq,  * hMaxbase;
  521.     char        ch;
  522.     long        len, maxlen, ibase;
  523.     short        spccount, nilcount;
  524.     short        maxlet, let, nseq, maxcnt;
  525.     short        letcount[kLetrange], letfirst[kLetrange];
  526.     float     rmaxc;
  527.  
  528.     DSequence* cons= this->Consensus();        
  529.     if (!cons) cons= (DSequence*) this->First(); //any seq will do...
  530.     
  531.     hMaxbase= NULL;
  532.     firstCommon= NULL;
  533.     if (cons) {        
  534.         maxlen= cons->LengthF();
  535.         hMaxbase= (char*) Nlm_MemGet(maxlen+1,true); 
  536.         firstCommon= (char*) Nlm_MemGet(maxlen+1,true); 
  537.         
  538.         for (ibase= 0; ibase<maxlen; ibase++) {
  539.             for (let= 0; let < kLetrange; let++) { letcount[let]= 0; letfirst[let]= kLetlast; }
  540.             nilcount= 0; 
  541.             spccount= 0;
  542.             nseq= 0;
  543.             
  544.             long iseq, maxseq= GetSize();
  545.             for (iseq=0; iseq<maxseq; iseq++) {
  546.                 DSequence* aSeq= this->SeqAt(iseq);
  547.                 if (!aSeq->IsConsensus() && aSeq->Kind() != DSequence::kOtherSeq) {
  548.                     hSeq= aSeq->Bases();
  549.                     len= aSeq->LengthF();
  550.                     if (ibase<len) ch= toupper( hSeq[ibase]);
  551.                     else ch= 0;
  552.                     if (ch >= 'A' && ch <= 'Z') {
  553.                         let = ch - 'A';
  554.                         letcount[let]++;
  555.                         if (letfirst[let] == kLetlast) letfirst[let]= Min(kLetlast, iseq);
  556.                         }
  557.                     else if (
  558.                          ch == DSequence::indelSoft 
  559.                         || ch == DSequence::indelHard 
  560.                         || ch == DSequence::indelEdge 
  561.                         || ch == ' ' || ch == '_' )  
  562.                         spccount++;
  563.                     else 
  564.                         nilcount++;
  565.                     nseq++;
  566.                     }
  567.                 }
  568.                 
  569.             for (let= 0, maxlet= 0, maxcnt= -1; let < kLetrange; let++) 
  570.                 if (letcount[let] > maxcnt) { 
  571.                     maxcnt= letcount[let]; 
  572.                     maxlet= let; 
  573.                     }
  574.                     
  575.             if (maxcnt>0) {
  576.                 firstCommon[ibase]= letfirst[maxlet];
  577.                 rmaxc= (100.0 * maxcnt) / nseq;
  578.                 if (rmaxc >= minCommonPerCent)   //mark maxlet's
  579.                     hMaxbase[ibase]= maxlet + 'A';
  580.                 else
  581.                     hMaxbase[ibase]= '!';
  582.                 }
  583.             else {
  584.                 firstCommon[ibase]= kLetlast;
  585.                 hMaxbase[ibase]= '!';
  586.                 }
  587.             }
  588.             
  589.         hMaxbase[maxlen]= 0;
  590.         firstCommon[maxlen]= 0;
  591.         }
  592.     
  593.     return hMaxbase;
  594. }
  595.  
  596.